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In the theory and practice of inverse problems for partial differential equations 

^^ (PDEs) much attention is paid to the problem of the identification of coeffi- 

cients from some additional information. This work deals with the problem 
of determining in a multidimensional parabolic equation the lower coefficient 
that depends on time only. To solve numerically a nonlinear inverse problem, 
^^ linearized approximations in time are constructed using standard finite element 

• procedures in space. The computational algorithm is based on a special decom- 

^ position, where the transition to a new time level is implemented via solving two 

'~~' standard elliptic problems. The numerical results presented here for a model 

I 2D problem demonstrate capabilities of the proposed computational algorithms 

^ for approximate solving inverse problems. 
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1. Introduction 

Mathematical modeling of many applied problems of science and engineering 
results in the numerical solution of inverse problems [2 [H [3] ■ Inverse problems 
often belong to the class of ill-posed (conditionally correct) problems, and there- 
fore various regularization algorithms are employed to solve them numerically 

SEE]. 

Particular attention should be given to inverse problems for PDEs [7l |8]. 
In this case, a theoretical study includes the fundamental questions of unique- 
ness of the solution and its stability both from the viewpoint of the theory of 
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differential equations [HI HO] and from the viewpoint of the theory of optimal 
control for distributed systems [TT]. Many inverse problems are formulated as 
non-classical problems for PDEs. To solve these problems approximately, em- 
phasis is on the development of stable computational algorithms that take into 
account peculiarities of inverse problems [HI [T^ . 

Among inverse problems for PDEs we distinguish coefficient inverse prob- 
lems, which are associated with the identification of coefficients and/or the 
right-hand side of an equation using some additional information. When consid- 
ering time-dependent problems, the identification of the coefficient dependences 
on space and on time is usually separated into individual problems [HI IH] ■ In 
some cases, we have linear inverse problems (e.g., identification problems for the 
right-hand side of an equation) ; this situation essentially simplify their study. 

Much attention is paid to the problem of determining the lower coefficient 
of a parabolic equation of second order, where, in particular, the coefficient 
depends on time only. An additional condition is most often formulated as a 
specification of the solution at an interior point or as the average value integrated 
over the whole domain. The existence and uniqueness of the solution of such an 
inverse problem and well-posedness of this problem in various functional classes 
are examined, for example, in the works [l4 | [T5 l [T6 t ll7j. 

Numerical methods for solving the problem of the identification of the lower 
coefficient of parabolic equations are considered in many works [T51 dni HOI HIl 
[21]. In view of the practical use, we highlight separately studies dealing with 
numerical solving inverse problems for multidimensional parabolic equations [231 
[24] . To construct computational algorithms for the identification of the lower 
coefficient of a parabolic equation, there is widely used the idea of transformation 
of the equation by introducing new unknowns that results in a linear inverse 
problem. 

In this paper, for a multidimensional parabolic equation, we consider the 
problem of determining the lower coefficient that depends on time only. Ap- 
proximation in space is performed using standard finite elements [25l [26] . The 
main features of the nonlinear inverse problem are taken into account via a 
proper choice of the linearized approximation in time. Linear problems at a 
particular time level are solved on the basis of a special decomposition into two 
standard elliptic problems. The paper is organized as follows. In Section [2| for 
a parabolic equation of second order, we formulate the inverse problem of the 
identification of the lower coefficient. The computational algorithm based on 
the linearization scheme is described in Section [3] Section [4| presents possibil- 
ities of the schemes with the second-order approximation in time. Numerical 
results for a model 2D inverse problem are discussed in Section [5] 

2. Problem formulation 

For simplicity, we restrict ourselves to a 2D problem. Generalization to the 
3D case is trivial. Let x — (a;i,a;2) and fi be a bounded polygon. The direct 
problem is formulated as follows. We search u{x,t), <t <T, T > such that 



it is the solution of the parabolic equation of second order: 

du 

— -dw{k{x)gia.du)+p{t)u = f{x,t), xefl, 0<t<T. (1) 

The boundary and initial conditions are also specified: 

k{x)--- + g{x)u = 0, x(^dn, 0<t<T, (2) 

u{x,0) — uo{x), X £ Q, (3) 

where n is the normal to fl. The formulation ([l])-(l3]) presents the direct problem, 
where the right-hand side, coefficients of the equation as well as the boundary 
and initial conditions are specified. 

Let us consider the inverse problem, where in equation (n]), the coefficient 
p{t) is unknown. An additional condition is often formulated as 

u{x,t)uj{x)dx^(p{t), 0<t<T, (4) 

Jn 

where uj{x) is a weight function. In particular, choosing uj{x) — d{x — x*) 
{x* G Q), where S{x) is the Dirac (5- function, from Q, we get 

u{x*,t) = ip{t), 0<t<T. (5) 

We assume that the above inverse problem of finding a pair of u{x,t), p(t) 
from equations ([l])-® and additional conditions Q or (l5]) is well-posed. The 
corresponding conditions for existence and uniqueness of the solution are avail- 
able in the above-mentioned works. In this paper, we consider only the nu- 
merical solution of these inverse problems omitting theoretical issues of the 
convergence of an approximate solution to the exact one. 

From the nonlinear inverse problem we can proceed to the linear one. Sup- 
pose 

v{x,t)=x{t)u{x,t), xW=exp( / p{0)de 
Then from ([l])-([3]), we get 

-^~div{k(x)gYadv)^x{t)f{x,t), xen, 0<t<T, 

dv 
k{x)— + g{x)v = 0, xedn, 0<t<T, 

v{x,0) — uo{x), X £ Q. 

The additional conditions Q and ((sl) to identify uniquely v{x,t), x(i) take the 
form 

v{x, t)uj{x)dx = xW'fiit), 0<t<T, 
n 



u{x*,t)^x{tMt), 0<t<T. 

The above transition from the nonhnear inverse problem to the Unear one is in 
common use for numerical solving problems of identification. In our work, we 
focus on the original formulation of the inverse problem Q-Q (or ([l])-(|3]), ([5|) 
without going to the linear problem. 



3. Computational algorithm 

The inverse problem of determining the pair of u{x,t), p{t) is nonlinear. 
The standard approach is based on the simplest approximations in time and 
involves the iterative solution of the corresponding nonlinear problem for the 
evaluation of the approximate solution at a new level. In our work, we apply 
such approximations in time that lead to linear problems for evaluating the 
solution at the new time level. 

Let us define a uniform grid in time 

w^ = w^ U {T} = {t" = Tir, n = 0, l,...,iV, tN^T} 

and denote y" = y{t"), i" = nr. Finite element approximations in space are 
employed. In the polygon i7, we perform a triangulation and introduce for this 
computational grid a finite-dimensional space V C H^{Q,) of finite elements. 

Using the fully implicit scheme for approximation in time, we obtain the 
following variational problem: 



/ vdx + / k{x)gradu"^^ giadvdx 

Jn ''' Jn 



+ / .g(a;)u"+iwd2:+p"+i / u''+\dx (6) 

Jon Jn 



J 

Jn 



f{x,r+')vdx, yveV, n = 0,l,...,7V-l, 



/ u^vdx = I uavdx. (7) 

Jn Jn 

The additional relations Q and (Is]) take the form 



u 



n+l, ,/'™\ J™ _ ,^"+1 



n 



uj{x)dx = (p"+\ (8) 



u 



n+l 



(a;*) = ^"+i, 7i = 0,l,...,iV-l. (9) 

To evaluate at the new time level the approximate solution u"+^(a;), p""*"^ 
from ([6])-([8]) or ^, ([t]), ^, some iterative procedures are necessary. In solving 
time-dependent problems, the solution slightly varies when it pass from the pre- 
vious time level to the next one. This basic feature of time-dependent problems 
is widely used in numerical solving nonlinear problems through the application 



of linearization procedures. We use a similar approach for the numerical solu- 
tion of the inverse problem that is concerned with the identification of the lower 
coefficient of a parabolic equation. 

Instead of (|6]), we will solve the following equation: 

/ vdx + / fc(a;) gradw""^^ gradwda; 

Jn '^ Jn 

+ [ g{xy+^vdx+p''+^ [ u"vdx (10) 

Jon Jn 

f{x,e+^)vdx, yveV, n = 0,l,...,N-l. 



In this case, the lower coefRcient of the parabolic equation is taken at the upper 
time level, whereas the approximate solution u{x,t) is treated at the previous 
time level. Let us consider the solution procedure of the problem ([7|, ([9]), (10 1 
in detail. 

For the approximate solution at the new time level m""*"^, we introduce the 
following decomposition [13l [27j : 



u"+i(a;) = y"+i(a;) + p^+^ w''+^ [x) . (11) 

To find y"+^(a;), we employ the equation 

vdx + I fc(a;) gradj/"^^ gradwda; + 

/on 

n = 0,l,...,iV- 1. 



/ vdx + / fc(a;) gradj/"^^ gradwda; + / g{x)y"'^^vdx 

Jn '^ Jn Jon 

= [ f{x,r+^)vdx, yveV, 
Jn 



(12) 



m 
The function w"^^{x) is determined from 

,n+l 



/ vdx + / /i;(a;) gradw"'*'"'^ gradrida; + / g{x)w^''^^vdx 

Jn ''' Jn Jan 

vJ'vdx, VweF, n = 0, 1,...,^-1. 



(13) 



n 



Using the decomposition (11|-(12), equation (10) holds automatically for any 

To evaluate p""^^, we apply the condition ([9]) (or ([8|). The substitution of 
( [TI| ) into ^ yields 



The fundamental point of applicability of this algorithm is associated with the 
condition w"+^(a;*) ^ 0. The auxiliary function iii"+^(a;) is determined from 



the grid elliptic equation ( 13 ). The property of having fixed sign for w"+^(x) is 



followed, in particular, from the same property of the solution at the previous 



time level u"(a;). Such constraints on the solution can be provided by the 
corresponding restrictions on the input data of the inverse problem. In any 
case, this problem requires special and careful consideration. In this paper, we 
assume that the constraint w"^^{x*) 7^ is satisfied. 



In solving problem (^8|-(10l, instead of (14), we have 



„n+l 



1 



,,ri+l 



uj{x)dx 



•^ 



n+1 



y'^+'uj{x)dx 



(15) 



under the condition that 



,,n+l 



ijj{x)dx 7^ 0. 



In this case, additional restrictions are formulated on the function a;(a;), e.g., 
its fixed sign in O. 

Thus, the computational algorithm for solving the inverse problem ([l| )-([4| 
(or ([l])-(||]), ([5|) based on the linearized scheme Q, (|8|, ([lO| (or (|7|,([9|, ( |lOf ) 
involves the solution of two standard grid elliptic equations for the auxiliary 
functions y"''^^{x) (equation (12)) and t(;"+^(a;) (equation (|13|)), the further 



the relation (11). 



evaluation of p"+^ from (15) (or (14|), and the final calculation M"+^(a;) from 



4. Scheme of the second-order accuracy 

The nonlinear inverse problem ([l])-(p]) is characterized by a quadratic non- 
linearity. When using the scheme with linearization ([7]), ([8]), (10), the nonlinear 
term is approximated with the first order with respect to t. It is possible to 
apply the linearized scheme of second order. Let us consider the approximation 

a(i"+i/2)^(^«+i/2^ ^ ^a(t"+i)6(t") + ia(t")6(t"+i) +0(r2). 

Approximation of equation ([I]) with the boundary conditions (pi) using the 
Crank-Nicolson scheme yields the linearized scheme 



,,n+l 



-vdx H — 



- / A:(a;) gradu"^^ graduda; H — / g{x)u^^'^^vdx 

2 Jo 2 Jgfi 

- / /i;(a;) gradu" gradwda; H — / g{x)u"vdx 
2 7n 2 Jg^ 



1 
2^ 



n+l 



u'^vdx + -p" 
n 2 



(16) 



,n+l 



vdx 



= f f{x,t"+^)vdx, \/veV, n = 0,l,...,N-l. 
Jn 



The scheme (l7|, ((sl), ( jlG] ) belongs to the class of linearized schemes. In com- 
parison with the scheme (l7|, ((sl), (10), it has a higher order of accuracy in 
time. 



To implement (I?]), (pi, (16), we again use the decomposition (11). In this 
case, for y"+^(a;), we have 

vdx + - / /c(a;) grady"^^ gradurfa; + - / g{x)y"^^vdx 

T 2 Jo 2 Jas;^ 



1 



A; (a;) grad u" grad vdx 



+ - g{x)u"vdx + -p" / y"+'^vdx 
2 7af7 2 Jfi 

= f f{x,t"+^)vdx, \fveV, n = 0,l,...,N-l. 
Jn 

The auxiliary function ^"^^(a;) is defined as the solution of the equation 

,,n+l 



(17) 



f w"~^ 1 f 

/ vdx ~\ — / /c(a;) grad w"^^ grad wda; 

Jn T 2 Jf2 



1 

2 Jdn 
1 



g{x)'w"^^vdx H — p 



,,n+l 



Tjda; 



(18) 



it"udx, WveV, n = 0, l,...,iV- 1. 



n 



Further, as in the case of the first-order scheme, we employ (15). 



The Crank-Nicholson scheme for numerical solving the direct problems for 
parabolic equations is not very often used in computational practice. It is 
inferior to the fully implicit scheme in sense of conservation of monotonicity 
(fulfilment of the maximum principle for the grid problem) , it has poor asymp- 
totic properties for solving problems with large integration time, and it is not 
unconditionally SM-stable scheme [551 [2S]- For this reason, it is appropriate 
to consider another variant of linearization of this inverse problem, where the 
second-order approximation is applied only for the nonlinear term. In this case. 



instead of (10), (16), we put 



vdx + / fc(a;) grad w"''" grad W(ia;+ / g{x)u^^ vdx 
n ^ Jn Jan 



+ ^P 



n+l 



u-'vdx + -p" / vJ^'^^vdx 



(19) 



f f{x,r+^)vdx, \fveV, n = 0,l,...,iV-l. 
Jn 



The numerical implementation of the scheme (ItI), (IS]), (19) is performed in the 
standard way using the decomposition (11). 



5. Numerical examples 

To demonstrate possibilities of the above linearization schemes for solving the 
problem of the identification of the lower coefficient of the parabolic equation. 



we consider a 2D model problem. In the examples below, we put 

k{x)^l, f{x,t)=0, uo{x)^l, xen, g{x)^lO, x e dn. 

The problem is considered on a triangular grid, which consists of 1,180 nodes 
(2,230 triangles) and is shown in FigjT] Here fi is the trapezoid with the vertices 
coordinates (0, 0), (0, 1), (1.5, 0.5), 1.5, 0). The calculations were carried out for 
T = 0.1. The coefficient p{t) is taken in the form 



p{t) = 



lOOOi, 
0, 



< i < iT, 
^T <t<T. 



(20) 



The solution of the direct problem (|T])-(|3| at the observation point is depicted 
in FigJ2] It was otained using the fully implicit scheme with different time steps. 
The solution at the final time moment is presented in FigjSJ 




Figure 1: The computational grid 



The results of solving the inverse problem with variuos grids in time are 
shown in Fig|4j The solution of the direct problem obtained with N = 1000 is 
used as the input data (the function ip{t) in the condition ^). It is easy to see 
that the approximate solution of the inverse problem converges with decreasing 
the time step. These results were obtained using the first-order scheme (It]), Q, 

Numerical results obtained for the above problem using the second-order 
scheme ([7|, Q, (16) are shown in FigJ5] For the discontinuous right-hand 
side (I20I), we observe characteristic wiggles of the identified coefficient. Such 



oscillations of the approximate solution are typical for the scheme (ItI), ^, (19). 

If the desired solution (the coefficient p{t)) is smooth, then the effect of using 

the second-order approximation is clearly expressed. As an example, we present 



— N = 25 

— N = 50 

— N = 100 

— N = 1000 




Figure 2: The solution of the direct problem at the point of observation 




Figure 3: The solution of the direct problem at i = T 




N = 25 
N = 50 
N = 100 
exact 



Figure 4: The solution of the inverse problem 




Figure 5: The solution of the inverse problem (the second-order scheme) 



10 



the results of numerical solving the inverse problem, where the lower coefficient 
(the exact solution) has the form 



P(i) = 



lOOOi 
1 + 500i2 ■ 



The approximate solution obtained via the first-order scheme (It]), (|9|, (fTol) is 
shown in FigM whereas FiglT] demonstrates the computations conducted by 
means of the second-order scheme ([t]), ([q]), (16 1. 



N = 25 
N = 50 
N = 100 
exact 




Figure 6: The scheme of first order (smooth solution) 
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